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I. INTRODUCTION 

Understanding the phase diagram of Quantum Chro- 
modynamics (QCD) is one of the major challenges of par- 
ticle physics. Running and upcoming, respectively, high- 
energy heavy-ion experiments at LHC, RHIC, SPS and 
FAIR, NICA and J-PARC call for a deeper understand- 
ing of hadronic matter under extreme conditions. The 
central issues are the equation of state of QCD and the 
nature of the phase transition from the confined hadronic 
phase with chiral symmetry spontaneously broken to the 
deconfined quark gluon plasma with chiral symmetry re- 
stored. The deconfinement phase transition is expected 
to be driven by the gluon dynamics, while the chiral 
phase transition, obviously, is due to strong interaction 
of the quarks, which, of course, is also mediated by the 
gluons. It is therefore by far non-trivial that both phase 
transitions occur at roughly the same temperature as ob- 
served on the lattice 

In quenched QCD the deconfinement phase transition 
is related to the center of the gauge group. Center sym- 
metry is realized in the low-temperature confined phase 
and spontaneously broken in the high-temperature de- 
confined phase . When quarks are included center sym- 
metry is explicitly broken and the deconfinement transi- 
tion is expected to become a cross-over. The features 
of the chiral phase transition dominantly depend on the 
quark masses, which explicitly break chiral symmetry, 
but also on the strength of the chiral anomaly. Further- 
more, lattice calculations also show that confinement is 
generated exclusively by the low-energy gluonic modes 
while also medium-energy gluon modes contribute to the 
order parameter of spontaneous breaking of chiral sym- 
metry, the quark condensate [3| . 

Meanwhile by various theoretical studies evidence has 
been accumulated that the phase structure of hadronic 
matter is much richer than originally thought [4] . As the 
chemical potential or baryon density is increased one ex- 



pects, on the basis of large- iVj, arguments, the existence of 
a "quarkyonic phase" where quarks and gluons are still 
confined but chiral symmetry is restored 0. At even 
larger chemical potential the transition to color super- 
conducting quark matter occurs with color-flavor locking 
Ifi;]. There are also speculations on the existence of a 
chiral critical endpoint for finite baryon density . 

Clearly the understanding of the phase structure and, 
in particular, the deconfinement phase transition re- 
quires, like the understanding of the confined phase itself, 
non-perturbative methods. A rigorous non-perturbative 
treatment of QCD is achieved on the lattice. The lattice 
method has been quite successful for quenched QCD but 
becomes extremely expensive when dynamical quarks are 
included and fails at large baryon densities due to the no- 
torious fermion sign problem at non- vanishing chemical 
potential. Therefore, alternative non-perturbative meth- 
ods which rely on the continuum formulation of QCD and 
thus do not suffer from the problem connected with the 
lattice treatment of fermions are very desirable. In recent 
years substantial progress has been made in first principle 
continuum QCD calculations, which rely on functional 
methods and do not suffer from the fermion problems of 
the lattice. Among these methods are Dyson-Schwinger 
equations (DSEs) in Landau (8|] and Coulomb gauge [9i], 
functional renormalization group (FRG) flow equations 
in Landau gauge [l^l and Hamiltonian Coulomb gauge 
[ill , and variational approaches to the Hamilton formu- 
lation of Yang-Mills theory in Coulomb gauge [12, [13 ■ 
These various continuum methods have all their own ad- 
vantages and drawbacks, and by combining them one 
expects to gain new insights into the non-perturbative 
regime of the theory. So far the deconflnement phase 
transition has been studied mainly in the FRG approach 
[11 in or in the DSE approach The deconfine- 

ment temperature was extracted either directly from the 
Polyakov loop |lH|, which (in the absence of quarks) is 
the order parameter of confinement, or from the dual con- 
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densate [l6|, which can be related to the Polyakov loop 

In this paper we study this transition in a variational 
approach to Hamiltonian Yang-Mills theory extending 
previous work 13] to finite temperature. 

In the variational approach at zero temperature the 
energy density is minimized using Gaussian t ype ansatze 
for the Yang-Mills vacuum wave functional jla [isj . So 
far, this approach has been mainly applied to study the 
infrared sector of Yang-Mills theory [13, [l^ , but was re- 
cently extended to full QCD In the present paper 
we extend this approach to finite temperatures and study 
the deconfinement phase transition. For this purpose we 
will minimize the free energy after making an appropri- 
ate ansatz for the density matrix. Some initial investiga- 
tions in this direction were undertaken in Ref. [211] , where, 
for simplicity, only the so-called subcritical solutions [l^ 
were considered, which do not produce confinement in 
the sense that they yield an infrared (IR) finite gluon en- 
ergy and not a linearly rising (static) quark potential. As 
a consequence the deconfinement phase transition could 
not be studied. Furthermore, in Ref. [2ll | the projection 
of the grand canonical gluon ensemble onto zero color 
was considered and it was found that the effect of color 
projection is negligible. We will therefore ignore color 
projection in the present paper. 

The organization of the paper is as follows: In Sec- 
tion |lT] we briefly summarize the basic ingredients of the 
Hamiltonian approach to Yang-Mills theory in Coulomb 
gauge, which are needed for its extension to finite temper- 
atures. In Section IIIII we introduce the grand canonical 
ensemble of Yang-Mills theory. We define here the com- 
plete basis of the gluon Fock space as well as the density 
matrix following Ref. [2l|. In Section ITVl we present the 
DSEs for the ghost and Coulomb propagator. In Sec- 
tion |V] we calculate the partition function, the entropy 
and the free energy. The variation of the free energy is 
carried out in Section IVII and the renormalization of the 
resulting DSEs is given in Section IVIII In Section IVIIII 
we present the results of an infrared analysis of the cou- 
pled DSEs and discuss the importance of the Coulomb 
term of the Yang-Mills Hamiltonian in Section IIXI Fi- 
nally, in Section IXII we present the numerical solution of 
these DSEs and determine, in particular, the deconfine- 
ment transition temperature. A short summary and our 

The details of the 



conclusions are given in Section IXIII 
IR analysis are presented in the Appendix. 



II. HAMILTONIAN APPROACH TO 
YANG-MILLS THEORY IN COULOMB GAUGE 

Below we briefiy summarize the basic ingredients of the 
Hamiltonian approach to Yang-Mills theory in Coulomb 
gau ge needed for its extension to finite temperatures 

To simplify the bookkeeping we will use the compact 
notation j4^J(xi) = A{1) for colored Lorentz vectors like 



the gauge potential and an analogous notation for col- 
ored Lorentz scalars like the ghost C°i(xi) = C(l). A 
repeated label means summation over the discrete color 
(and Lorentz) indices along with integration over the d 
spatial coordinates 

A-B^A{1)B{1)^ / d'xJ2 E A"(x)i?r(x). (1) 

•' i=l a=l 

We define in coordinate space 

6{l,2)^6^'^H,,,,{xi)S{^,-X2), (2) 

where 

did. 



(3) 



is the transverse projector. Furthermore, indices will be 
suppressed when they can be easily restored from the 
context. 

After resolving Gauss's law in Coulomb gauge the 
Yang-Mills Hamiltonian 0| reads in our notation 



H=\[JAMl)JAm) + B{l)B{l)\ 
+ Y^aV(1)^-4^a(1,2)p(2) 



(4) 



= Hk + He 



He, 



where H(l) = —i5/5A{l) is the canonical momentum 
(electric field) operator, B{1) is the non-abelian magnetic 
field 

B° = V X A" + I /"''^A'' X A" (rf = 3), (5) 
g is the coupling constant, and 

G:,1(1,2) = H'^^'^^a^^ -5Arnxi)ar)'5(xi -X2) 

(6) 

is the Faddeev-Popov determinant with A"'' = jacb^c 
being the gauge field in the adjoint representation of the 
gauge group SU(iVc) with structure constants Fur- 
thermore, 



p(l) = p"i(xi) = i?(l;2,3)A(2)H(3), 



(7) 



i?(l; 2, 3) = r^"-^"^ <5,,,3 (5(xi - X2) <5(xi - xg) 

is the color charge density of the gluons and 

Fa{1, 2) EE F^i"^(xi, X2) = Ga{1, 3) Go\3, 4) G^(4, 2) 

(8) 

is the so-called Coulomb kernel, with Gg ^ being the bare 
inverse ghost operator, obtained by setting ^ = in 
Eq. ^ . In the presence of matter fields with color charge 
density Pm(l), the gluon charge p(l) in the Coulomb term 
He [Eq. dH)] is replaced by the total charge p(l) + Pm(l) 



3 



and the vacuum expectation value of ^"^(1,2) acquires 
the meaning of the static non-abehan Coulomb potential. 

The gauge fixed Hamiltonian Eq. (|4|) is highly non- 
local due to Coulomb kernel i^y!i(l,2), Eq. ([S]), and due 
to the presence of the Faddeev-Popov determinant J^, 
Eq. In addition, the latter also occurs in the func- 
tional integration measure of matrix elements of opera- 
tors O between Coulomb gauge wave functionals 



(V'i|0|^2)= I DAJAri[A]OMA.] 



(9) 



where the integration runs over transverse field config- 
urations. In practical calculations the elimination of 
unphysical degrees of freedom via gauge fixing is usu- 
ally beneficial (and sometimes unavoidable), in spite of 
the increased complexity of the gauge-fixed Hamiltonian 
Eq. . The non-trivial Faddeev-Popov determinant re- 
flects the intrinsically non-linear structure of the space of 
gauge orbits and dominates the IR behavior of the the- 
ory. Once Coulomb gauge is implemented, any functional 
of the (transverse) gauge field is a physical state. 



III. THE GRAND CANONICAL ENSEMBLE OF 
YANG-MILLS THEORY 

We are interested in the behavior of Yang-Mills theory 
at finite temperatures. For this purpose we consider the 
grand canonical ensemble of Yang-Mills theory, which 
is defined in the Hamiltonian approach by the density 
matrix 



V = exp(-^iJ) 



(10) 



with /3 — 1/fcgT being the inverse temperature mea- 
sured in units of the Boltzmann constant fcg . Formally, 
V [Eq. (fTO)) ] looks like the canonical ensemble since the 
chemical potential of the gluon vanishes. However, the 
thermal expectation values 



TiVO _ ^fc(fc|P0|fc) 



(11) 



have to be taken over the whole Fock space, i.e. over 
states {|fc)} with an arbitrary number of gluons. Since 
H is non-linear and non-local it is clear from the very 
beginning that we have to resort to approximations. We 
will use an analogous approximation scheme as in the 
zero-temperature case, see Ref. 13]. 

The trace in Eq. ([TT|) can, in principle, be calculated 
in any complete basis {|fc)}. Since we are mainly inter- 
ested in {{H)) we will choose a basis which is adapted to 
the structure of the Yang-Mills Hamiltonian. Following 
Refs. [l3| and [2l| we choose the basis of the gluonic Fock 
space of the form 



{A\k) = J^'/'{A\~k), 



(12) 



where \k) denotes a complete set of states of the Yang- 
Mills Fock space to be specified later. This ansatz re- 
moves the Faddeev-Popov determinant from the inte- 
gration measure in Eq. In this basis the thermal 
expectation value Eq. pT|) reads 



^^^^^ ^ E.mm 



(13) 



where we have introduced the abbreviation 

O = j\^^OJ^^^^, (14) 
which yields for the density matrix Eq. ([TO)) 



(15) 



The transformed Yang-Mills Hamiltonian H defined by 
Eq. is obtained by replacing in the Yang-Mills 

Hamiltonian Eq. (|4]) the momentum operators J^^'^HJ^H 
by Htfi, whereEl 



ii = jl/W/^ = n + A^ 

^ ^ 2 SA 



(16) 



To work out the thermal averages it is convenient to go 
to momentum space 



A{x) = J dfce"'-^A(k), 
H(x) = J dfce-*-^H(k), 



(17) 



(dfc = d'^fc/(27r)'^) where the Fourier components satisfy 
the commutation relation 



[A(1),H(2)] =15(1,2). 



(18) 



In the standard fashion we express the Fourier compo- 
nents in terms of creation and annihilation operators 



^i(k) 



i=(<(k)+af(-k)), (19a) 



n?(k) = i^^ (af (k) - <(-k)) , (19b) 

which satisfy the usual commutation relations 

[a^(k),a^^^(q)] =5'^''i,,(k)(2^)'^5(k-q), (20) 



^ Note that since the color density of the gluon field p JT]! is linear 
in the momentum operator 11, the Coulomb Hamiltonian |(4]| 
has, concerning the occurrence of EI and J a, the same structure 
as the kinetic term ~ J^^JIJ^II. This is not surprising since 
the Coulomb term is the longitudinal part of the kinetic energy 
^ J^^nll Jylllll with n" being the longitudinal momentum op- 
erator determined by the resolution of Gauss's law. 



4 



where iij(k) 



k/|k| is the transverse 



projector Eq. in momentum space. Here w(k) is, so 
far, an arbitrary (positive definite) function of momenta. 
It enters the vacuum state |0) defined by 

<(k)|6)=0, (21) 

which has the "coordinate" representation 

{A\6) = AAcxp [-iA(l)c.(l, 2)A{2)] , (22) 

where is a normahzation constant. The states 



|0), afikm, ar(k)ar(q)|0). 



(23) 



form a complete basis in the gluonic Fock space for any 
positive definite function w(k). Below we will use this 
basis to evaluate the thermal expectation values Eq. (IT51) . 

Because of the non-local structure of the Yang-Mills 
Hamiltonian [Eq. (jj])], the density matrix Eq. (fTS)) can 
only be treated in an approximate fashion. To make the 
actual calculation feasible we follow Ref. [2l| and replace 
the full Yang-Mills Hamiltonian H in the density opera- 
tor 2? [Eq. dTSl) ] by a single-particle operator 



h = 



(24) 



where J7°''(k) will be later determined by minimizing the 
energy density, yielding 



V = cxp(-/3/i). 



(25) 



By global color and rotational invariance, 57"'' (k) is color- 
diagonal and independent of color, f2°''(k) = S°'^il{k), 
and furthermore depends only on |k|. The same is true 
for the kernel w in the vacuum wave functional Eq. (|22p . 
Since the transformed density matrix 2? [Eq. (|25p] is the 
exponential of a single-particle operator Wick's theorem 
applies to the thermal averages Eq. ([T^. whose tem- 
perature dependence is exclusively given by the finite- 
temperature Bose occupation numbers n{k), defined (in 
momentum space) by 



(a^t(k)a5(q)) = S^%{k){27rfS{k~q)n{k), 
n{k) ^ [exp(/317(k)) - 



(26) 



From Eqs. (|19ap and ()19bp one finds for the gluon prop- 
agator 

Z?(1,2):=(A(1)A(2)), i^(k) = i±M^, (27) 



2w(k) 



for the momentum propagator 



if (1,2) (H(1)H(2)), Kik) = i±|^a;(k). 



and for the "mixed" propagator 

(A(1)H(2))^^<5(1,2), 
the latter being independent of the temperature. 



(28) 
(29) 



IV. GHOST AND COULOMB PROPAGATOR 

We will resort to the same approximation used in 
Ref. in the T = case, i.e. calculating the energy 
up to two loops. To this order the following representa- 
tion of the Faddeev-Popov determinant holds j24| 



where 



JA = exp[-A(l)x(l,2)A(2)], 

X(l 2) = -i/^!l^\ 
' ' 2 \ 5A{1)6A{2) I 



(30) 



(31) 



is the so-called curvature, which, in fact, represents the 
ghost loop (see below). Strictly speaking this represen- 
tation was derived only for T = 0. However, the proof 
given in Ref. [IJ] can be straightforwardly extended to 
finite temperatures provided the density matrix T) is the 
exponential of a single particle operator, so that Wick's 
theorem holds for the thermal averages (O) [Eq. (IT3|) ]. 
This is the case for the density matrix Eq. With 
the inverse Faddeev-Popov operator Ga [see Eq. dH])] the 
curvature Eq. (pij) can be expressed as 

X(l, 2) - i (ro(l; 3, 4) G^(3', 3) ro(2; 4', 3') Ga(4, 4') 



where 



ro(i;2,3) 



M(l) 



(32) 



(33) 



is the bare ghost-gluon- vertex. Defining the ghost prop- 
agator by 

G = (Ga) (34) 
and the full ghost-gluon vertex F by 

{GaToGa) = GFG, (35) 
the ghost loop Eq. (15^ becomes 

X(l, 2) = iro(l; 3, 4) G(3', 3) r(2; 4', 3') G(4, 4'). (36) 

Furthermore, the ghost propagator G satisfies the Dyson 
equation [l^j 

G = Go + GoSG, (37) 

where 

Gf{^,y) = [i-d'rX'^ (38) 

is the free ghost propagator and 

S(l, 2) f o(3; 1, 4) G(4, 4') D{3, 3') f (3'; 4', 2) (39) 

is the ghost self-energy. As shown in Landau gauge [1^ 
the ghost-gluon vertex is not renormalized; this applies 
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also in Coulomb gauge, see Ref. [l8|. At zero temperature 
the full ghost-gluon vertex F can be, to good approxi- 
mation, replaced by the bare one Fq [1^. We assume 
that this approximation works also at finite tempera- 
tures. Assuming a bare ghost-gluon vertex and express- 
ing the ghost propagator Eq. p4p in momentum space 
by the ghost form factor 



(40) 



the curvature Eq. p6p becomes in momentum space 

and the ghost form factor (i(k) satisfies the following DSE 

1 



where 



/d(k) =iV,y"dq[l- (k.q)2] 



d(k - q) H- 2n(q) 
(k-q)2 2c.(q) 



(42) 



(43) 



As at zero temperature we express the Coulomb prop- 
agator 



i^(l,2) = (i^^(l,2)) 



(44) 



by means of the Coulomb form factor /(k) defined in 
momentum space by 

9^F^\k) = gH^^ G(k)/(k) k2G(k) = 5-^^ /(k). 

k 

(45) 



Using the operator identity 



(46) 



and, proceeding as in the zero-temperature case (see e.g. 
Ref. [13]), neglecting the (/-dependence of the density ma- 
trix one finds the relation 



,2 9 



/(k) = V^rf-'(k) 



(47) 



from which one finds with Eq. (I42p the following integral 
equation 



/(k) = l + /^-(k) 



(48) 



W^N.j^,[l-ii..c,)Y-^^^^^^ (49) 

1 + 2n(q) 
2c^(q) 

This equation is formally the same as at zero temperature 
except that the gluon propagator is replaced by its finite- 
temperature counterpart, Eq. (j27p . 



V. THE FREE ENERGY 

At vanishing chemical potential the thermodynamic 
potential of the grand canonical ensemble is given by the 
free energy 

T={H) - TS, (50) 
where S is the entropy, which is defined by 



S — — fcg Tr 



V , V 
— m — 

Z Z 



where 



Z = TrV 



(51) 



(52) 



is the partition function of the grand canonical ensemble 
with V defined by Eq. (fTS)) . By straightforward algebraic 
manipulation the entropy Eq. (I5ip can be converted to 



S — \nZ 



/3 



dlnZ 
dp 



(53) 



In principle, we could calculate the free energy from the 
partition function, Z = exp(— /3J^). For the exact den- 
sity matrix Eq. (jlSp this would yield the same result as 
Eq. (j50p . However, in the present case, where we have 
replaced the full density matrix Eq. (fTS)) by the one of 
a system of independent quasi-particles, Eq. (j25l) . it is 
mandatory to evaluate J- from Eq. (|50p in order to cap- 
ture the essential correlations between the gluons. Be- 
sides, the density matrix T) Eq. ([25)) has so far not been 
determined. 

The partition function of the density matrix 2? 
[Eq. (^5]) ] is calculated in the standard fashion yielding 



Z = exp 



(d - l){Nl - l)V / dk ln(l -f n(k)) 



(54) 



where d — \ — taili) is the number of independent po- 
larization directions in d spatial dimensions and N"^ — 1 
is the number of color degrees of freedom of the gauge 
bosons. With Eq. ([M)) we find from Eq. (l55|) for the 
entropy density s[n\ per degree of freedom of the gauge 
bosons, defined by 



1 



S={d-l){N',-l)Vj—^s[n], (55) 



the expression 



s[n] ^ j dfc [/3n(k)r2(k) + ln(l + n(k))] , (56) 

where 57 (k) is the variational kernel in the density ma- 
trix [Eq. (gH)]. To obtain the free energy from Eq. (j50p 
we have still to calculate the energy (if), which is done 
below. 
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In the evaluation of the expectation value of the Yang- 
Mills Hamiltonian Eq. (0]) we use the same approximation 
as in the T = case keeping only terms up to two loops in 
the energy [Tsl . The magnetic energy is straightforwardly 
evaluated using Wick's theorem. One finds 



V / dp-^[l + 27i(p)] 
g'N^iN^ - 1)^ I ^^^^ d(d - 3) + 3 - (p . q)2 



16 



^(p)'^(q) 

[1 + 2n(p)] [1 + 2n(q)] . (57) 



The evaluation of the kinetic energy and the Coulomb 
energy is somewhat more involved. With the representa- 
tion Eq. we obtain from Eq. (fTOj) 



n(i) = n(i)-ix(i,2)A(2). 



(58) 



For the expectation value of the two momentum opera- 
tors entering the kinetic and Coulomb terms one finds by 
using Eqs. ^ and (P7 |) - (pg)) 

(flt(l)n(2)) = K{1, 2) 2) l')D{l', 2')x(2', 2). 

(59) 

Contracting the external indices one obtains from this ex- 
pression immediately the expectation value of the kinetic 
term 

= i [K{1, 1) - 1) + X(l, 2)i?(2, 3)x(3, 1)] 

+ x'(q)^(q)]- (60) 

This expression is not obviously positive definite (as it 
should be); however, it is not difficult to show that this is 
indeed the case. Using Eqs. (P7| and (^5]) and separating 
the zero- and finite-temperature terms, Eq. (1601) can be 
written as 



Hk) = 



+ n{q) 



2c^(q) 

^'(q) + x'(q) " 

w(q) 



(61) 



from which it is seen that (Hk) is indeed positive definite. 

Restricting ourselves to including up to two (over- 
lapping) loops in the energy allows us to replace the 
Coulomb kernel Fa (1,2) in He by its expectation value 



Eq. (j44p . resulting in the approximation 
(p(l)F^(l,2)p(2)) 

~ {FAih 2))i?(l; 3, 4)i?(2; 5, 6)(A(3)nt(4)A(5)n(6)) 
= F(l, 2)i?^(l; 3, 4)i?,(2; 5, 6) [(^(3)^(5)) (fit (4)n(6)) 

+ (A(3)n(6))(nt(4)A(5))]. (62) 

The remaining expectation values can be straightfor- 
wardly carried out using Eqs. ([55)) and ([5^ . Then we 
find for the Coulomb term 



(He 



|-F(l,2)i?(l;3,4)i?(2;5,6) 



[oiS, 5) [if (4, 6) - x(4, 6) + x(4, 4')i?(4', 6')x(6', 6)] 
J(3,6)-i?(3,3')x(3',6) 
"i^(4,5)-x(4,4')i?(4',5)]}. (63) 



As noticed in Ref. [21], with the replacement 
<(x,y)^<5'^'' 
the Coulomb Hamiltonian Hq [Eq. ^] becomes 



He 



(64) 



(65) 



where 



(66) 
(67) 



is the total color charge and we have used 

[Q", Ja] =0. 

The last relation holds since the Faddeev-Popov determi- 
nant Ja is invariant under global color rotations, which 
are generated by Q". This is also explicitly seen by using 
the representation Eq. ([30p . In a colorless universe 

((Q^g'^)) = (Q'^Q'^) = iQ^Q'^) ^ (68) 

holds. To ensure that this condition ()68l) is respected by 

(He) in Ref. HH the Coulomb kernel ^Ff{x,y) was 
replaced by 

^2 
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y^^A (x,y) 



(69) 



Furthermore, in Ref. [21[ the Yang-Mills grand canoni- 
cal ensemble was projected onto zero total color. It was 
found in there that the effects of both color projection 
and of the replacement Eq. (j69p is negligible. Therefore 
in the following we will ignore the projection on zero color 
states [as well as the replacement Eq. ([M]) ]. 

Rewriting the above given thermal average of the 
Hamiltonian, Eqs. ([57)) . (150)) and (1551) . in momentum 
space one finds for the energy density per degree of free- 
dom e[7i,cj] defined by 

V , , 



{H) = {d-l){N',-l) 



{2nY ^ ' J' (70) 



the following expressions 



ec 
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,K=^^ldq [Kid) - x(q) + (q)^(q)] , (71a) 
es^^^ /d..^i.(.) + ^!«^ /d.d,Z.(pM,), (71b) 
ec = ^0^(2^)' J dpdq[d-2 + ip- q)2]F(p + q){i^(p)if (q) + i?(q)i^(p) + i?(p)i?(q) [x(p) - x(q)] ' 

+ [Z?(p)-Z?(q)][x(p)-x(q)]-i}, (71c) 



where i^(k) is the Coulomb propagator defined by 
Eqs. (gH) and (gS]). 



VI. THE FINITE-TEMPERATURE 
VARIATIONAL PRINCIPLE 

The kernel r2(k) defining the density matrix "D [see 
Eqs. and ^S^l] is so far completely arbitrary. We 
will now determine it by the finite-temperature varia- 
tional principle. At fixed temperature and volume and 
arbitrary particle number the thermodynamic potential 
to be minimized is the grand canonical potential, which 
in the present case coincides with the free energy since 
the chemical potential of the gluons vanishes. Therefore 
we minimize the free energy or its density 



f[n,uj] = e[n,u] - Ts[n], 



(72) 



where e[n,CLi] is the energy density [Eq. (j70p ] and s[n\ is 
the entropy density [Eq. ([55]) ]. Instead of varying the free 
energy density with respect to r2(k) it is more convenient 
to take the variation with respect to the occupation num- 
bers n(k) [Eq. which is equivalent since n(k) is a 
monotonic function of f2(k). From Eq. (|56p we find for 
the variation of the entropy 



(73) 



Therefore stationarity of the free energy density, 5f/5n = 
0, requires 



r2(k) 



5e[n^ uj] 
(5n(k) 



(74) 



which, in the spirit of Landau's Fermi liquid theory, iden- 
tifies i7(k) as the quasi-gluon energy. Of course, this re- 
sult could have been anticipated from the form of the 
finite-temperature Bose occupation numbers, Eq. (|26p . 

So far the kernel w(k), which defines the vacuum wave 
functional [Eq. (P^ ] and thus our basis of the Fock space, 
is completely arbitrary and we could use any positive defi- 
nite kernel and the corresponding gluon basis to calculate 
the thermodynamic averages. As long as we include the 
complete set of states and keep the full canonical density 



operator the thcrmodynamical averages are independent 
of i^(k). Thus, in principle, the free energy should not 
depend on w(k). However, due to the truncation of the 
full Hamiltonian in the density operator T) [Eq. ([T5|)] to 
a single particle one [Eq. ([25]) ] the actual choice of the 
basis, i.e. of w(k), does matter and the optimal choice of 
w(k) is obtained by extremizing the free energy 



(5a; (k) 



0. 



(75) 



For the evaluation of 5e[n, a;]/(5a;(k) we skip the implicit 
w(k) dependence of x(k) and F{k), since their inclusion 
would give rise to higher-order loops. Ignoring this im- 
plicit a;(k) dependence the energy density depends on 
n(k) and a;(k) only through the gluon field and momen- 
tum propagators, D [Eq. ^] and K [Eq. (EH])], i.e. 



j] = e [Z3[7i,a;], i^[n, w]] 



(76) 



Using the chain rule and the explicit form of the propa- 
gators [Eqs. (HZl) and ^Ml] we obtain 



5e 



1 



5e 



(5n(k) 

5e 
(5a; (k) 



a;(k) (5L>(k) 
1 + 2n(k) 



+ a;(k) 

5e 
SD{k) ' 



5e 



SK{k) ' 
-^2(k) 



de 



6K{k) 



(77) 
(78) 



Since the entropy density s[n] [Eq. (1551) ] does not explic- 
itly depend on aj(k) the condition Eq. ((75|) reduces to 

6e[n,io]/SujCk) = 0, (79) 

which in view of Eq. (j78l) leads to the condition 

Se n,, , Se 



dD{k) 



SK{k) 



(80) 



For T = 0, which implies n(k) = 0, this equation reduces 
to the gap equation obtained in Ref. 13]. Inserting the 
gap equation dHHl) into Eq. ([771) we find from Eq. 



n(k) = 2a;(k) 



6e 



SK{k) 



With the explicit expressions for the energy density given 
in Eq. ([TT]) we obtain 



n{k) 



l + /o(k)^ 



(82) 
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where 



(83) 

Using the expHcit expressions for the energy density 
[Eq. (ITT]) ], the gap equation can finally be expressed as 



^2(k)=k2 + x2(k)+/0[n]+/^[n](k), 



(84) 



where 



= (85) 



2d 



is the tadpole term stemming from the non-abelian part 
of the magnetic energy, and 

^^N(k)^^/d.[.-2 + (k.,)^]M 

X {(1 + 2n(q)) [..^(q) _ ^2^y^ ^ (^^^^ _ ^(k))'] 

-2c^(q)(x(q)-x(k))} (86) 

is the contribution of the Coulomb Hamiltonian. These 
loop integrals are ultraviolet (UV) divergent and re- 
quire regularization and eventually renormalization of 
the gap equation. Fortunately the temperature depen- 
dence of these loop integrals (which is due to the finite- 
temperature occupation numbers n{\i.)) does not give 
rise to additional UV singularities. Therefore the zero- 
temperature counterterms are, in principle, sufficient to 
eliminate the UV singularities. However, special care is 
required to separate the temperature dependence from 
the UV-singular terms, which will be done in the next 
section. 



VII. RENORMALIZATION 

At large momenta |k| ^ k^T the temperature should 
become irrelevant. Consequently, the finite-temperature 
solutions w(k), d(k), x(k) should possess the same UV 
behavior as in the zero-temperature case. Indeed, the 
finite-temperature contributions to the loop integrals 
(the terms proportional to the occupation number n(k)) 
are all ultraviolet finite. This is because for |k| — )■ oo we 
have r2(k) ~ w(k) ~ k and thus the finite-temperature 
occupation numbers 



n(q 



oo 



-/3n(q— >oo) 



-I3q 



(87) 



cut off the large momenta. Therefore, the renormaliza- 
tion can be done independent of the temperature sub- 
tracting only zero-temperature counterterms. 

For the renormalization we follow Ref. [Tsj and sepa- 
rate the various degrees of UV divergences of the loop 
integrals of the gap equation by writing Eq. (j86p as 

/^[n](k) = /(f)(k,r) + 2x(k)/W(k,T)-l-/„M(k), (88) 



where the integrals 



2{d-l) 



d(7[d-2 + (k-q)2] 



i^(k~q) 
w(q) 



Mq)-X(q))'-Hk)-X(k))'l (89) 



are linearly (for 1 = 1) and quadratically (for / ~ 2) UV 
divergent, while the finite-temperature contribution 



J.[n](k) = 0^ j dg[d-2 + (k.q)2] 



j"(k - q) 
w(q) 



X n(q) uj\q) - ^^(k) + (^(q) - x(k))' (90) 



is UV convergent due to Eq. (l87l) . Analogously we write 
for the tadpole Eq. ([55]) 



(91) 



The loop integrals I^^ (k, T) and /° are defined as in the 
zero-temperature case, Ref. [l^. However, their entries 
a;(k), d(k), and x(k) are temperature dependent. These 
integrals contain all UV divergences, while the finite- 
temperature modifications /i^[n](k) and /^[n] are UV fi- 
nite. Therefore, the renormalization can be carried out, 
in principle, in the same way as in the zero-temperature 
case, Ref. [27]. However, due to the implicit temperature 
dependence of w(k), (i(k), x(k) we have to be careful not 
to introduce finite-temperature effects by the renormal- 
ization. Following the renormalization procedure given in 
Ref. [131 for T = and subtracting only zero-temperature 
terms one arrives at the renormalized gap equation 



c^2(k) = k^ + x2(k) + A/(f)(k) + A/^ + /0[r 



Co 



2x(k)[A/«(k) + ci] +/.N(k), (92) 



where 



X(k) = x(k) - x(Mx)It=o 



A/W(k) = 4')(k) 



(0) 



A/" 



(93) 



and Co, c\ are finite renormalization constants surviving 
from energy counterterms [l^ 

Ai/ = icoA(l)A(l) + iCiA(l)H(l). (94) 

Co and c\ are the finite parts of the divergent constants 
Co and Ci, respectively, i.e. Ci = Cf™ + Ci. Furthermore, 
/ijj, is an arbitrary scale arising from the renormalization 
of the Faddeev-Popov determinant 22|. The subscript 
T = in Eq. (j93p means that the corresponding quan- 
tities have to be taken with the zero-tempcrature self- 
consistent solution. The choice of these finite renormal- 
ization parameters will be discussed in Sect. [X] 
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The renormahzed ghost DSE is obtained from Eq. 
with a subtraction at an arbitrary scale fid &t T = 

d-i(k) = rf-i(/i,)-A/rf(k), 

AIdik)^Id{k)-Id{fid)\j.^o- 
The Gribov-Zwanziger confinement scenario requires the 
horizon condition 



d~\k = 0) 



T=0 



= 0, 



(96) 



to be satisfied at T = 0. This condition can be explicitly 
built in the renormahzed ghost DSE (j95l) by choosing the 
renormalization constant d{iJ,d) such that 

d^^ifid) = A/d(fc = 0)\r^^Q = Id{0)\j,^Q ~ /d(Aid)|y^o , 

(97) 

which can be fulfilled for arbitrary renormalization scale 
fid- Inserting this value into Eq. (j95p we obtain 

d-i(k) = -(/rf(k)-/,(0)|^^J, (98) 

which is nothing but the renormahzed ghost DSE ([55)) 
with the renormalization scale fixed at /i^ = [and the 
horizon condition Eq. (j96p built in] . This shows that im- 
plementing the horizon condition automatically puts the 
renormalization scale in Eq. ([M]) to /x^ = 0. Solutions of 
the coupled DSEs and gap equation satisfying the hori- 
zon condition are called "critical" and those with 
d-^{0) > subcritical [22|. 

The renormahzed DSE for the Coulomb form factor 
/(k) [Eq. is analogously obtained by subtracting 

Eq. (j48l) once at T = and an arbitrary renormalization 
scale fif, yielding 



/(k) 
A//(k) 



(99) 



/(/i/) + A//(k), 

^/(k)-//(M/)|T=o- 
In principle, the renormalization scale /i/ of the Coulomb 
form factor can be independently chosen from that of 
the ghost fid- However, since /(k) and d{'k) are tightly 
related by Eq. (|47| for consistency one should choose 
fif = f^d- 

Within our approach the finite-temperature Yang- 
Mills theory is now determined by the following set of 
coupled equations: The Eq. ([5^ for the quasi-gluon en- 
ergy f](k), the gap equation ^ for w(k), the DSE ^ 
for the ghost form factor d(k) and the DSE (IMl) for the 
Coulomb form factor /(k). 

A comment is here in order: In principle, we would 
get only two equations from the finite-temperature varia- 
tional principle, one equation for 57 (k) and one for w(k). 
However, in the process of evaluating {H) we have in- 
troduced, for convenience, additional propagators (ghost 
and Coulomb propagator), which we did not explicitly 
express as functionals of uj{k) (and i7(k)). Instead of 
that we derived DSEs for these quantities, which all con- 
tain loop integrals. In taking the variation of the energy 
{H) with respect to w(k) and 51 (k) the implicit w(k) — 
and J7(k) — dependence of these propagators is ignored 
since it would give rise to two-loop terms in the equations 
of motion, see Ref. [23 for more details. 



VIII. INFRARED ANALYSIS 

Before presenting the numerical solutions of the cou- 
pled Eqs. dsn, dHll), dMl) and dHU), we investigate their 
infrared behavior. At zero ternperature the infrared anal- 
ysis was carried out in Refs. O, [11]. At arbitrary finite 
temperature the infrared analysis of the DSEs cannot 
be done in the usual way. This is because the finite- 
temperature occupation numbers n(k) [Eq. (|26l) ] depend 
exponentially on the gluon quasi-energy f2(k). The in- 
frared analysis can, however, be carried out in the high- 
temperature limit where the occupation numbers become 



n(k) 



ML 
fifk) 



(100) 



This limit is sufficient to exhibit the infrared behavior of 
the propagators in the deconfined phase. With the rep- 
resentation Eq. (jlOOp the infrared analysis of the coupled 
Dyson-Schwinger equations can essentially be carried out 
as at T = 0. 

In general, the IR analysis can be carried out in two 
ways: i) using the angular approximation [l^ replacing 
kernels iir(k— q) depending on the difference between the 
external momentum k and the loop momentum q by 

if (k - q) ^ K{k) e(fc -q)+ K{q) Q{q - k), (101) 

with k — |k|, q = |q|. Although being approxima- 
tive, this method has the advantage that IR power law 
ansatze for the propagators are indeed restricted to the 
IR regime, ii) Without using the angular approxima- 
tion the IR power law ansatze have to be used in the 
loop integrals for the whole momentum range [l^, [29l |. 
This method is, in principle, exact in the IR (up to the 
omission of possible logarithms) but has the disadvan- 
tage that the UV behavior of these integrals is usually 
changed by the IR power law ansatze and UV-finite loop 
integrals may turn into UV-divergent ones. Fortunately 
both methods yield very similar results as we will explic- 
itly demonstrate in the Appendix. 

For the infrared analysis we assume the following 
power law ansatze 



w(k ^ 0) 



A 
A;" 



d(k ^ 0) 



B 



(102) 



For r = the infrared analysis with and without angu- 
lar approximation was carried out in Refs. jl3|] and [l8|, 
respectively. Assuming the horizon condition [Eq. ([96])] 
and a bare ghost-gluon vertex one finds from the ghost 
DSE (I42p in d spatial dimensions in both cases the sum 
rule 



a = 2/3 + 2 - d. 



(103) 



Including also the gap equation the infrared exponents 
are fixed and one finds, abandoning the angular approx- 
imation, in d = 3 spatial dimensions two solutions |18(] 



/3 ~ O.i 



/3 



1. 



(104) 
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(In the angular approximation only the second solution is 
obtained.) Both exponents are also found in the numer- 
ical solutions, and were originally obtained in Ref. [l^l 
(/3 ~ 0.8) and in Ref. (/3 = !)• In d = 2 spatial di- 
mensions one finds a single solution with /3 = 1/2 in the 
angular approximation and (3 = 0.4 when the angular ap- 
OTOximation is abandoned, while the numerical solution 
1^] yields /3 ~ 0.44. 

For high temperatures using the approximation 
Eq. (jlOOp the infrared analysis is carried out in the Ap- 
pendix. From the ghost DSE one finds the same sum rule 
Eq. (|103p as in the zero-temperature case. Including also 
the gap equation one obtains in d = 3 at high tempera- 
tures a solution with the infrared exponent of the ghost 
form factor 

/3 = 1/2 (105) 

both with and without the angular approximation. With 
this value for /3 one finds from the infrared analysis of 
the DSE for the Coulomb form factor /(k -> 0) 
the infrared exponent Awl which together with /3 = 
1/2 for the ghost form factor leads to a linearly rising 
Coulomb potential F(k — >• 0) ~ l/fc'', as is also found 
in the lattice calculation [3lj |. The infrared exponents 
obtained in the infrared analysis are confirmed by the 
numerical calculations, as we will see in Section EH 

IX. NEGLECT OF THE COULOMB TERM 

We are mainly interested in the description of the de- 
confinement phase transition. If the Gribov confinement 
scenario is realized in Coulomb gauge, this transition 
should manifest itself in a change of the infrared behav- 
ior of the various propagators, in particular in that of 
the ghost form factor d(k) and the gluon energy w(k) 
and ri(k), respectively. In the (renormalized) gap equa- 
tion ([^^ the contributions from the Coulomb term are 
infrared sub-leading ( A/(') (fc = 0) = 0) even for infrared 
finite a;(k). Furthermore, the ghost DSE does not re- 
ceive contributions from the Coulomb term. Therefore, 
for a first qualitative description of the deconfinement 
phase transition the Coulomb term should be negligible. 
With the neglect of the Coulomb term, the gap equation 
reduces to 

(k) = k^ + x' (k) + A/^ + 7^ [n] +c, + 2cix(k) , (106) 

while the ghost DSE ((98|) remains unchanged. At T = 
both AI^ [Eq. dM])] and I^[n] [Eq. dH])] vanish and the 
gap equation p06p reduces to 

c^2(k)=.k2 + x2(k) + 2cix(k)+co. (107) 

Figure [1] shows the result of the numerical solution of the 
ghost DSE and the full gap equatior0 (|92]) at T = and 



^ When the Coulomb term is included the DSE for the Coulomb 
form factor has also to be solved. This was done as described in 



the corresponding solutions of the gap equation (|107p 
with the Coulomb term neglected for Cq = ci = 0. 
As is seen, there are only small deviations in the mid- 
momentum regime, while both solutions agree in the ul- 
traviolet and, in particular, in the infrared. Furthermore, 
in a quasi-particle description of the gluon sector (under- 
lying the present approach) the neglect of the Coulomb 
term is conceptually advantageous (and expected to give 
a better description) as long as two particle correlations 
are neglected. The reason is the following: The Coulomb 
term gives rise to the UV-singular quasi-gluon self-energy, 
Ar2(k) = w(k)/n(k), see Eq. ([5^ . which is diagram- 
matically illustrated in Fig. [5^. In a two-(quasi-)gluon 
state these divergent self-energy contributions are pre- 
cisely canceled by the divergent contribution from the 
Coulomb interaction to the two-gluon energy, shown in 
Fig. Therefore it does not make sense to keep the 
Coulomb term in the quasi-gluon energy as long as the 
two-body correlations are not taken into account We 
will therefore neglect the Coulomb term in the following. 
Then the loop integral in Eq. ((5^ has to be discarded so 
that 

f](k)=w(k). (108) 

Furthermore, A/° [Eq. ([M]) ] represents the change of 
the tadpole due to the change of the self-consistent so- 
lution a;(k) at finite temperature (relative to the zero- 
temperature case) but does not contain the change of 
the tadpole due to the explicit temperature dependence 
of the gluon propagator via the finite-temperature occu- 
pation numbers. Therefore, we expect AI^ to be small 
and we will neglect it. Then the finite-temperature gap 
equation p06p reduces to 

Lo^ (k) = k^ + (k) + /o [n] + 2ci x(k) + cq . (109) 

Here only /^[n] (the temperature-dependent part of the 
tadpole, see Eq. ((9T|) ') is explicitly temperature depen- 
dent, while the curvature x(k) depends only implicitly on 
the temperature via the ghost form factor. However, as 
we will see in Sect. IXI[ it is this temperature dependence 
of the ghost which triggers the deconfinement phase tran- 
sition. 



X. CHOICE OF RENORMALIZATION 
CONSTANTS 

Let us now discuss the choice of the finite renormaliza- 
tion parameters. Since the renormalization can be com- 
pletely accomplished by renormalizing the theory at zero 



Ref. [l3ll replacing d(k) by its bare value <i(k) = 1 in the loop 
integral //(k). 

^ Note also that in QED the Coulomb term vanishes identically in 
the absence of external charges. 
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FIG. 2. Contributions of the Coulomb term to (a) the gluon 
self-energy and (b) the two-gluon energy. A wavy line rep- 
resents the gluon propagator, the double line stands for the 
Coulomb kernel. 



temperature, the remaining finite-temperature renormal- 
ization constants cq , ci should be also fixed at T = and 
then kept fixed independent of the temperature. To fix 
these constants we notice that the T = solutions are 
IR divergent with uj{k) ^ xi^) ~ ^/k", a = 0.8, 1, see 
Sect. I Villi For such solutions the gap equation (|107p re- 
duces in the IR with Eq. (fTU^ to 

L.(fc)-x(fc) = ci + ||A;" + ... (110) 

The constant c\ obviously determines the infrared limit 
of w(k) — x(k), and ci = is required in order that the 
't Hooft loop obeys an area law, Ref. [131 . This value 
is also favored by the variational principle (of minimal 
energy), Ref. [l^]. Furthermore, recent lattice investiga- 
tions of the Yang-Mills vacuum wave functional in c? = 2 
spatial dimensions show that for this value of c\ the wave 
functional Eq. (j22p yields statistical weights for abelian 
plane-wave configurations which are in very good agree- 
ment with the ones of the exact vacuum wave functional 
[33 |. We will therefore put ci = 0. 

As is seen from Eq. (|110l) the renormalization constant 
Co is IR subleading compared to c\ and influences the 
mid-momentum regime, which is however expected to 
affect the deconfinement phase transition. Fortunately 
in the present variational approach cq needs not to be 
treated as a free parameter but can be determined by 
minimizing the energy density. First notice that with 
ci = the IR behavior of the T = gap equation (|110p 



reduces to 

^(fc)-x(fc) = ||fc" + --- (111) 

Furthermore, with Eq. (I30p . which is correct to the two- 
loop order (in the energy) considered in the present pa- 
per, our vacuum wave functional defined by Eqs. \Y1\ 
and (|22p reads 

(A|0) - AAexp ^(1) H\ 2) - x(l, 2)]^(2)| . 

(112) 

In order that this wave functional is regular | (AjO) p < oo 
for all A we must have 

w(fc)-x(fc)>0 (113) 

According to Eq. (|llip this requires for small fc 

Co > 0. (114) 

We will see that this condition is also required by the 
positivity of the energy. Furthermore the numerical re- 
sults given in Sect. IXII will show that the energy density 
takes its minimal value for co = 0. 

Using the gap equation (jl07l) the zero-temperature 
energy density Eq. (j71j) with the counterterm ~ Co 
[Eq. (IMl) ] fully includecO can be cast into the form 

e = eK + eB = ^ydA;[c^(k)-x(k)]. (115) 

Note that the energy density does not explicitly but 
merely implicitly depend on co, since w(k) and x(k) de- 
pend on Co via the gap equation (|107p . From Eq. (|115p 
it is seen that stability of our non-perturbative vacuum 
requires again the condition (jllSp . i.e. cq > 0. 



* The counterterm ~ Ci [Eq. 1941 1] is only needed for the Coulomb 
energy. 
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FIG. 3. Energy density [Eq. (|115|l ] as function of the renor- 
malization constant Co for the solution with IR exponent 
a = l. 



Figure [3] shows the energy density e [Eq. (|115p ] as 
function of cq for the solution with IR exponent a = 1. 
The energy density Eq. (jllSp takes its minimal value at 
Co = 0, which is therefore the physical value and thus 
chosen in alia variational calculations. 

There is one hidden renormalization constant left, 
which is the scale fj,^ in the subtracted curvature xO^) 
Eq. (j93|) . This quantity can be chosen in a wide range 
without significantly changing the resulting propagators, 
see Ref. [1J|. (Clearly, this parameter cannot be chosen 
fiy- = since this would result in w(0) — 0, which is in 
contradiction to the Gribov scenario and also with the 
findings on the lattice (SSj).) 



XI. NUMERICAL RESULTS 

The coupled ghost DSE (|M| and the gap equation 
(I109p can be solved numerically for the whole momentum 
range by iteration in the way described in Ref. Il£| for the 
zero-temperature case. (For details see also Ref. [ll|.) In 
the following we will restrict ourselves to d = 3 spatial 
dimensions and to Nc = 2. For the numerical calcu- 
lation it is convenient to introduce dimensionless quan- 
tities, rescaling all dimensionful quantities with appro- 
priate powers of an arbitrary momentum scale /i. The 
rescaled dimensionless quantities will be indicated by a 
bar: 



k:=- 



aJ(fc) 



etc. 



(116) 



Because of the use of a logarithmic momentum grid an 
infrared cut-off Air, is needed, which we choose, in dime- 
nionless units, in the range Air — 1Q~^ . . . 10~^. 



temperature case, from which the physical scale is fixed. 
We keep the remaining undetermined renormalization 
constant fixed a,t jl^ = 4, but its precise value is irrelevant 
for the qualitative behavior of the obtained solutions0 

The Gribov-Zwanziger confinement scenario requires 
d^^{k = 0) = 0. Solutions which satisfy this condi- 
tion are referred to as "critical" , while solutions with 
d~^{k = 0) > are called "subcritical" . Figure H] 
shows the results of the numerical solutions of the ghost 
DSE dnSI and the gap equation pIB)) at T = and 
various choices of the ghost renormalization constant 
d~^{Hd = Air), starting at d~^{Xi^) > 3Air and succes- 
sively decreasing d~^{XiR). As long as £i~^(AiR) ^ 2Air 
the solutions obtained are subcritical, i.e. d~^{0) > 0. 
At sufficiently small d~^(AiR) w 2Air the critical solu- 
tion with /3 ~ 0.8 is obtained. When c?~^(Air) is further 
decreased critical solutions with (3 — 1 appears. 

The critical solutions can be more easily found by im- 
plementing the horizon condition d~^{k = 0) = explic- 
itly into the ghost DSE (IMI) by choosing its renormaliza- 
tion constant d~^{fid = 0) = [resulting in Eq. ([55]) ] 
and solving the coupled equations and (llOOp it- 

eratively. Choosing as starting point of the iteration 
(^"^(Air) ^ Air one obtains the critical solution with 
j3 ~ 0.8, while the choice d~^(AiR) ^ Air results in a 
continuous set of critical solutions with (3 = 1 differing in 
the value of the IR coefficient B of d(k), see Eq. (|102p . 
Let us also mention that the convergence of the itera- 
tion is generally faster when the initial functions d(k) 
and a;(k) satisfy the IR sum rule Eq. (jl03l) and the re- 
lation (jA.37p . which represents the IR limit of the ghost 
DSE, see the Appendix. For /3 = 1 and rf = 3 Eq. ((X37)) 
reduces to 



A 
52 



87r2 



(117) 



In Fig. [S] the infrared coefficients B [Eq. (I102p ] of the 
various /3 = 1 solutions have been rescaled to make them 
coincide in the IR, so they differ in the mid-momentum 
regime. 

As already mentioned before our equations can be en- 
tirely expressed in terms of dimensionless quantities, re- 
flecting the scale invariance of Yang-Mills theory. To 
confront our numerical results with the experimental or 
lattice data we have to fix the (so far arbitrary) scale /i, 
see Eq. (|116p . For this purpose we use the lattice data 
for the gluon propagator, which were obtained by fix- 
ing the scale by means of the Wilsonian string tension 
CTw — (440MeV)2. (This quantity has not yet been cal- 
culated in the present Hamiltonian approach.) The lat- 
tice calculations carried out in Coulomb gauge in Ref. [1^ 



A. Zero-temperature solutions 

Before presenting the finite-temperature solutions, it 
is worthwhile and necessary to reconsider the zero- 



^ Note that, in principle, the value of this finite renormalization 
constant could also be determined by minimizing the free energy 
with respect to this parameter. This is, however, numerically 
extremely expensive and beyond the scope of the present paper. 
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U\i k/H 
FIG. 4. The ghost form factor (left panel) and the gluon energy (right panel) for different values of d~^(AiH,) at T = 
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FIG. 5. The different critical solutions with /3 = 1 overlayed to the same infrared coefficient B. 



show that the gluon energy aj(k) can be nicely fitted by 
Gribov's formula (sj) 

uj{k) = ^k' + -^, (118) 

with M ~ 880 MeV ~ 2y^. The Gribov formula im- 
plies the IR exponent a — 1 [cf. Eq. (I102p ] and hence 
by the IR sum rule (|103p /3 = 1. Therefore, to fix our 
scale from the lattice result for M, we have to use the 
/3 = 1 solutions. As discussed above, there is a whole set 
of /3 = 1 critical solutions differing in the IR coefficient 
of the ghost form factor B = B/^. From these we chose 
the one which fits Gribov's formula Eq. (jll8p best. For 
the above adopted values JI^ — 4 and cq = one finds 
for this solution the value 

S = — = 6.01 ±0.05. (119) 
M 

From the Gribov formula Eq. (I118P one reads off the IR 
coefficient (jl02p of the gluon energy 

A = M'^. (120) 

With this result we find from Eq. (jll7p for our scale 



This fixes our physical scale in terms of the Gribov mass 
M. 



B. Finite-temperature solutions 

In this section we use the same procedure described 
above to solve the coupled equations and (|109p for 
finite temperatures. The renormalization described in 
Sect. IVlIl reauires subtractions with the zero-temperature 
solutions in order to have temperature independent 
renormalization constants. However, such subtractions 
are numerically unstable. Because of this, in the nu- 
merical results presented in this section the subtrac- 
tion of the ghost loop has been done with the finite- 
temperature solutions. Consequently, the horizon con- 
dition (i^^(/c = 0) = is built in explicitly at any tem- 
perature. 

Starting with the critical solutions at T = and in- 
creasing the temperature the self-consistent solutions re- 
main more or less unchanged up to a critical tempera- 
ture Tc, where the IR exponent of the ghost form factor 
suddenly changes to /? « 0.5 and the gluon energy be- 
comes infrared finite [a ~ 0) in accord with the sum rule 
Eq. (jl03p . This is nicely seen in Fig. [S] where we show the 
infrared exponent for the critical solutions as function of 
the temperature. The two critical solutions with /? ~ 0.8 
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FIG. 6. Infrared exponent of the ghost form factor as function 
of the temperature for the critical solutions at T = with [3 ~ 
0.8 (full line) and /? = 1 (crosses). Above Tc both solutions 
merge to a single solution. 
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FIG. 7. Solutions of the ghost DSE for different temperatures 
T. 



and P — \ existing below Tc merge to a single solution 
which approaches /3 = 0.5 for T S> Tc. The IR exponent 
P — 0.5 is precisely the value of the ghost form factor in 
d — 2. Thus, the change of the IR exponent at Tc is in 
accord with dimensional reduction. 

Figures [7] and [5] show the self-consistent finite- 
temperature solutions for the ghost form factor and the 
gluon energy as function of the momentum for various 
temperatures. As the temperature is increased above 
Tc the plateau value of the gluon energy starts increas- 
ing linearly with the temperature as is seen in Fig. [HI 
where we show the infrared value a; (Air) as function of 
the temperature. The linear increase with the tempera- 
ture of a;(AiR) is easily understood by noticing that above 
Tc the temperature is the only energy scale and for di- 
mensional reason a;(AiR) has therefore to scale with the 
temperature. 

The infrared value of the gluon energy w(Air) can be 
interpreted as an effective gluon mass. Zooming into the 
behavior of a;(AiR) near Tc, which is done in Fig. [TOl 
we can extract the critical exponent k of this quantity 
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FIG. 8. Solutions of the gluon gap equation for different tem- 
peratures T. 




T/T, 



FIG. 9. tj(AiR) as function of the temperature. 



defined by 

From our numerical solution we extract a value of k 
0.37, which is similar to the value of k ~ 0.41 obtained 
in Ref. 35], where a quasi-gluon picture has been used 
to fit the lattice results for the energy density and the 
pressure, and furthermore using critical exponents from 
the d — 3 Ising model, which is in the same universality 
class as SU(2) gauge theory. 

The sudden change of the ghost infrared exponent, see 
Fig. [5] (or alternatively the drop in the IR value of the 
gluon energy ^(Air), see Fig. ^ can be used to deter- 
mine the critical temperature Tc- From our numerical 
solutions we extract with Tc ~ fiTc and Eq. (I12ip for 
SU(2) the value of 

275---290MeV, (123) 

which is somewhat smaller than the lattice result of Tc ~ 
295 MeV [SU(2)]. 

The value of the critical temperature is seen to be in- 
sensitive of the actual value of < d~^{ijLd = Air) ~ Air. 
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FIG. 10. Critical behavior of the eflective gluon mass for 
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FIG. 11. aj(AiR) as function of the temperature for various 
values of the couphng g. 

Yet there is one parameter which can effect the form of 
the phase transition: the couphng constant g entering the 
tadpole term I^in]. This term acts as an effective mass 
in the gap equation (I109|) . As long as the temperature is 
small, we expect the IR behavior of the gluon energy to 
be insensitive of the precise value of this effective mass. 
However, when the temperature is raised the effect of the 
mass scale cannot be neglected. When we increase the 
value of g (Fig. [TT|) the phase transition is weakened. For 
the self-consistent solutions we have chosen g — 0.5. 

Removing the tadpole, the free QED solution (i^(k) — 
k, a = —1 and P — 0) exists at high temperatures. 

In Fig. 1121 we show the running coupling constant a{k) 
calculated from the ghost-gluon vertex, as described in 
Ref. fis"], for both T = and T > and normalized to 
the infrared value ac (at T — 0). In both cases a{k) is 
IR finite but the IR plateau value decreases by an order 
of magnitude above the deconfinement phase transition. 
Below and above Tc this plateau value remains more or 
less unchanged as the temperature varies. The IR finitc- 
ness of the running coupling is guaranteed by the sum 
rule Eq. (|103p of the IR exponents of ghost and gluon 
propagators, which holds both at T = and T — > oo. 
As is well known, in the high-temperature limit d = 4 




0.0001 0.001 0.01 0.1 1 10 100 1000 



FIG. 12. Zero- and finite-temperature running coupling con- 
stant normalized to the infrared value ac at T = 0. 



100 



10 - 



1 r 




0.00010.001 0.01 0.1 1 10 100 1000 



FIG. 13. Gluon propagator D(k) [Eq. ((27|)] at zero and finite 
temperature. 



Yang-Mills theory reduces to o? = 3 Yang-Mills theory 
coupled to a Higgs field, which is the temporal compo- 
nent of the gauge field in d = 4. Since the Higgs field 
certainly contributes to the running coupling constant 
there is no obvious reason why at high temperatures the 
running coupling should approach that of the (confining) 
d = 3 Yang-Mills theory, although it is IR finite in both 
d = 4 Yang-Mills theory at T — > oo and d = 3 Yang- 
MiUs theory at T = 0. 

Finally, in Fig. [13] we show the gluon propagator 
Eq. (P7)) at zero and at finite temperatures above the 
deconfinement phase transition. Below the deconfine- 
ment transition temperature the gluon propagator is in 
accord with Gribov's formula (jll8l) vanishing in the IR, 
while above the deconfinement phase transition the gluon 
is massive. Furthermore the effective IR gluon masses 
increase with the temperature as is explicitly seen in 
Fig. [11] Note also that in the UV the gluon propaga- 
tor is independent of the temperature. 
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XII. SUMMARY AND CONCLUSIONS 

In this paper we have studied the grand canonical 
ensemble of Yang-Mills theory in the Hamiltonian ap- 
proach in Coulomb gauge and investigated the finite- 
temperature deconfinement phase transition. For the 
density operator a quasi-particle picture was assumed 
and the quasi-particle energies were determined by min- 
imizing the free energy. A complete basis of the glu- 
onic Fock space was constructed by creating an arbi- 
trary number of quasi-gluons on top of a Gaussian-type 
vacuum wave functional whose width was considered as 
variational kernel and determined by minimizing the free 
energy. This results in the finite-temperature gap equa- 
tion, which has to be solved together with the Dyson- 
Schwinger equations for the ghost and Coulomb propa- 
gators. We have shown that the effect of the Coulomb 
term of the Yang-Mills Hamiltonian (which represents 
the longitudinal part of the kinetic energy) is negligible. 
Neglecting this term the gap equation and DSE for the 
ghost propagator decouple from the DSE for the Coulomb 
propagator. We have solved these equations analytically 
in the high-temperature limit in the infrared regime. We 
have found that the infrared exponents of the gluon en- 
ergy (a) and of the ghost form factor {(3) satisfy in the 
high-temperature limit the same sum rule a = 2/3 ~ 1 
as at zero temperature. While at zero temperature two 
solutions with /3 = 1 and (3 ~ 0.8 exist, in the high- 
temperature limit there is only a single solution /3 = 1/2. 
Our numerical solution of the coupled ghost DSE and gap 
equation shows that at a critical temperature Tc there is a 
sudden change of the infrared exponents from their zero- 
temperature values to their high-temperature limits. At 
this deconfinement phase transition both T — critical 
solutions with (3 ~ 0.8 and (3 = 1 merge to a single so- 
lution with /3 = 1/2. In accord with the sum rule at the 
critical temperature the gluon energy changes from being 
infrared divergent {a ~ 0.6 and a = 1) to being infrared 
finite (a = 0). This shows that while the gluons are ab- 
sent from the infrared spectrum in the confined phase 
they become massive particles in the deconfined phase. 
For the effective gluon mass we have found a critical ex- 
ponent of 0.37. From the sudden change of the infrared 
exponents and of the infrared value of the gluon mass we 
have extracted a critical temperature of the deconfine- 
ment phase transition of Tc = 275 . . . 290 MeV using the 
lattice results for the Gribov mass. 

An alternative way to describe the deconfinement 
phase transition is to study the behavior of the Polyakov 
loop , as was done in the ERG approach in Landau gauge 
UM- This can be also done in the present approach and 
will be subject to future work. 

Towards the extension of the present approach to full 
QCD one has first to study the gauge group SU (3), which 
has a first order phase transition. From the lattice stud- 
ies carried out in Ref. [s^ in Landau gauge one may 
expect that the order of the phase transition manifests 
itself in the critical behavior of the effective gluon mass 



close to the phase transition, but this is an open issue 
and requires further studies. 

The results obtained in the present paper are rather 
encouraging for an extension of the present approach to 
full QCD at finite temperature and baryon density. A 
first step in this direction was recently undertaken by 
studying chiral symmetry breaking in a variational ap- 
proach to QCD in Coulomb gauge at zero temperature 
and density |20] . The extension of this approach to finite 
temperature and non- vanishing chemical potential is, in 
principle, straightforward but technically involved. Since 
this approach is formulated in the continuum it will not 
face the problems one encounters in the lattice formula- 
tion at finite chemical potential. 
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Appendix: Infrared analysis 

Below we carry out the IR analysis of the coupled DSEs 
in the high-temperature limit. As usual we assume power 
law ansatze in the infrared (fc = |k|) 

ABC 
'^ik) = —, d{k) = -^, ^^^^^^^ 

and, in addition, for the ghost form factor d{k) the hori- 
zon condition 

d-\k = Q)=Q, (A.2) 

so that (3 > Q. The coefficients A, B. . .as well as the 
IR exponents a, /3. . . are expected to depend on the tem- 
perature. The modifications at finite temperature arise 
exclusively from the extra piece of the gluon propaga- 
tor Eq. ([?7)) containing the finite-temperature occupa- 
tion numbers n{k) Eq. (|26p . defined in terms of the ker- 
nel ri(/c) occurring in the ansatz for the density matrix, 
Eq. Obviously, for £ > the finite-temperature 

part of the gluon propagator is IR subleading since in 
this case 

n{k ^ 0) = (A.3) 

and we expect the zero-temperature result for the sum 
rules of the IR exponent to remain true at finite T. Note 
also that the property (jA.3|) is maintained in the high- 
temperature approximation Eq. (|100l) . 
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1. Infrared analysis in angular approximation 

We begin with the analysis of the Dyson-Schwinger 
equation for the ghost form factor. It differs from the 
zero-temperature equation only by the replacement of 
the zero-temperature gluon propagator l/[2a;(fc)] by its 
finite-temperature counterpart [l-t-2n(fc)]/[2w(fc)]. As in 
the r = case it is convenient to consider the derivative 
of the ghost DSE ^ 



dk 



d-\k) = i'M- 



(A.4) 



In the angular approximation, Eq. (jlOip . after the angu- 
lar integration over the d — 1-sphere 



1 



d'^-if7[l-(k-q)2] = 



1 



1 



(47r)'i/2 r(l + I) 



(A.5) 



we obtain for the loop integral 

Nr. d-l 1 (dikW 



I'M = 



2 (47r)'i/2 r(l f ) V / 



R{k), (A.6) 



where 



R{k) 



dqq'^-'^Diq) 



k 1 

d-l 



dqq" 



uj{q) 



(A.7) 



and we have used here the high-temperature limit 
Eq. (llOOp . For small k we can safely use the infrared 
asymptotic forms (jA.ip for a;(fc) and Vl{k) in the inte- 
grand, which yields 



R{k) 



A{a + d) [ E 



a 



(A., 



Inserting this result into (|A.6p we find from (IA.4I) the 
relation 



A _Nr d 



1 



1 



/3 + 2 



2 {A^TY/^T{l + !^)P{a + d) 



1 



E 



^Q-2/3+d-2 

a + d 



a + e d 



(A.9) 



where the left-hand side is a constant for fixed T. For 
e > the second term in the bracket is not IR leading 
compared to the first one and we obtain the relation 



2/3 - d -I- 2 , 



(A.IO) 



which is just the zero-temperature infrared sum rule [ISj. 
In the opposite case, e < the second term is IR leading 
and the sum rule 



e^2f3-d + 2 



(A.ll) 



follows. Note the two sum rules merge continuously at 
e = 0. 



Consider now Eq. (15^ for the gluon quasi-energy 
Obviously, if the Coulomb term is neglected, we have 
il(fc) — uj{k) and consequently e — a. For the infrared 
analysis of the full equation it is again convenient to take 
the derivative 



dfc \u{k)) 

With the ansatze (jA.l[) the left-hand side yields 
dfc \uj{k) J ' 



(A.12) 



(A.13) 



Using the angular approximation (jlOip for the Coulomb 
kernel F(k — q) we find for the derivative of the loop 
integral 



g^Nc 



2(47r)d/2 r(l + I) 



F'{k)R{k) 



(A.14) 



where R{k) has been defined by Eq. (|A.7I) . Using the 
explicit form of the Coulomb kernel Eq. pS)) 



g'F{k) 



d{kff{k) 



(A.15) 



we obtain with the ansatze (lA.ll) the infrared behavior 



F{k) 



(A.16) 



With this expression and (IA.8I) we find from Eq. (|A.14p 

g^Nc d - 1 B^L 



Inik) 



2(47r)'^/2r(l + f) A 
2/3 + 2 + X 



^a-2^+d-3-A 



a + d 



1 + 



k^T a + d 
'e~ a + e + d 



k' 



(A.17) 



Inserting Eqs. (|A.13p and (|A.17p into (IA.12P we find for 
£ > the relation 



while for e < 



e = 2(3-d + 2 + X, 



2e = 2l3'd+2 + X 



(A.18) 



(A.19) 



holds. Using the results (lA.lOp and (jA.lip from the ghost 
equation we find in both cases the sum rule 



a + X. 



(A.20) 



In the case of an IR finite Coulomb form factor (A — 0) 
the latter equations implies 

£ = a, (A.21) 

so that il(fc) and a;(fc) have the same IR behavior, up to 
possible logs. Indeed with A = and e = a > satisfying 
the sum rule (jA.lOp we find from Eq. (jA.17p 



Inik) 



(A.22) 
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implying that n{k) has an extra infrared logarithm in 
addition to the infrared power law inherited from uj{k) [f] 



n{k) =w(fc)(l+-lnfc), 



fc ^ 0. 



(A.23) 



The expression (j4ip for the curvature x{k) is formally the 
same as in the zero-temperature case and consequently 
its infrared analysis can be carried out as in the zero- 
temperature case, Ref. fl3| . In d spatial dimensions one 
obtains from the derivative of Eq. (PT|) for the curvature 
with the ansatze (IA.1|) 



C 



N, 



1 



/3- 



52 2 (47r)rf/2 r(l + d/2) 7(d - ;3) 
from which we find 

7 = 2;3 - d + 2 , 



^d_2/3-2+7^ (A.24) 



(A.25) 



which is the zero-temperature result. Inserting this rela- 
tion into eqs. (lA.lOp and (jA.llI) we obtain 



a = 7, 
a + e — J, 



e > 
e < 0. 



(A.26) 



Consider now the finite-temperature gap equation 
Proceeding as in the zero-temperature case [ij] one easily 
shows that for IR divergent x(fc), i.e. 7 > 0, the gap 
equation (jA.ip reduces in the IR limit to 

a;2(fc->0) =X^(fc^O) (A.27) 

implying that in the ansatze (jA.ip we have 

A^C. (A.28) 



With this equality the left-hand sides of Eq. (jA.9p 
(ghost DSE) and of Eq. (IA.24P (curvature) become equal. 
Equating the right-hand side of these equations and us- 
ing the sum rules ()A.10p and (IA.25P to cancel the IR 
power-laws we obtain the relation 



/3 



/3- 



la(d-/?) P{a + 2) 



(A.29) 



Using the sum rule (jA.lOp again to eliminate a in favor 
of (3 we finally obtain 



1 = 



+ 1) 



d- 1 {2(3 + 2 - d){d - P) ■ 
The solutions of this equation are 

fi^l id = 3), 



(A.30) 



(A.31) 



Note that the Coulomb integral lQ{k) [Eq. Il83l l] is positive def- 
inite since its integrand is positive, so that f2(fc), like u]{k), is 
positive definite. 



which were quoted in Sect. I Villi If the curvature x(^) 
is IR finite, i.e. 7 < 0, the IR dominant terms on the 
right-hand side of the gap equation (|92|) is non-zero due 
to the temperature-dependent part of the tadpole, I^[n]. 
Therefore in this case uj{k) is IR finite and non- vanishing, 
implying a ~ 0. With this value the sum rule (jA.lOp 
yields 



d-2 



I.e. 



/? = - (d = 3). 



(A.32) 



(A.33) 



This is the solution found in the numerical calculations 
in the deconfined regime, while the solutions (|A.3ip are 
realized in the confined phase. 



2. Infrared power law ansatze in loop integrals 



Using the power law ansatze (jA.ip for the whole mo- 
mentum regime the loop integrals (^5]) . ([M)) and (|106p 
can be calculated analytically. For the zero-temperature 
case this has been done in Refs. [l^ and 22 1. With 
the high-temperature limit of the occupation numbers 
Eq. (jlOOp these results can be extended to finite temper- 
atures. In this subsection, for simplicity, we will neglect 
the Coulomb term entirely, so that il{k) — uj{k). 

Using the infrared ansatze (jA.ll) one finds from the 
ghost DSE gl]) 



A_ ^ iVc(d-l) 
52 4 



r,d-2-|-Q-2/3 



(47r)''/2 



2k T 

hu> (a, /3, d) -t- fc" -^Idu {2a, (3, d) 



{AM) 



where 

Iduj{a,l3,d) 



r(i + f)r(i-f)r(i(2-d-a + /3)) 
r(i-f)r(d + f-f)r(i + |) 



(A.35) 

The new (finite-temperature) element is the second term 
in the bracket of Eq. (r04| . At T = or for a > (and 
arbitrary T) this term is IR subleading and one finds the 
sum rule 



a = 2/3-d + 2 



(A.36) 



already previously obtained in angular approximation, 
see Eq. (lA.lOp . With this sum rule (and a > 0) 
Eq. (|A.34I) reduces in the IR to 



A 
52 



Nc{d- 1) 



1 



IdLj{a,l3,d) . 



4 (477)^^/2 

Putting here /3 — I and d — 3 we obtain 

A_^N^ 
52 87r2 ' 



(A.37) 



(A.38) 
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which is the relation (|117p quoted in Sect. IXI Al 

For T > and a < the second term in Eq. (|A.3ip 
becomes IR dominant and we find the sum rule 



2a = 2(3 - d + 2 , 



(A.39) 



which was also obtained in the previous subsection in the 
angular approximation. 

We now turn to the gap equation (jl06p . The neglect of 
the Coulomb terms is irrelevant for its IR analysis since 
these terms are IR subleading as discussed in the previous 
subsection. 

Using the power law ansatz (|A.ip for the ghost form 
factor in the whole momentum regime one finds for the 
ghost loop Eq. (gl]) 



4 (47r)'^/2 



Idd{(3,d), (A.40) 



where 



IddiP, d) 



r(|-|m/3 + i-|) 
r(i + §)2r(d-/3) 



(A.41) 



From (jA.40l) we read off the IR exponent Eq. (lA.ll) of 
x{k) to be given by 



7 = 2^ - d + 2 , 



(A.42) 



which is again the same result as obtained in the angular 
approximation (|A.25p . We can distinguish now two cases: 

i) If the ghost loop x(^) is IR divergent, i.e. 7 > 0, 
the gap equation (11061) reduces in the IR to 



(A.43) 



and we find 7 = 0;, which also follows by combining 
(jA.36|) and (jA.42|) and agrees again with the results 
obtained in the angular approximation. 

Since a;(fc), x(fc) > 0, Eq. (IA.43P reduces then to 



(A.44) 



Inserting here for xi^) the expression (|A.40p and 
using the sum rule (|A.36p one finds the relation 



A _ Nc{d-2) 
52 ^ 8 



(47r)^/- 



:Iddil3,d). 



(A.45) 



Combining this equation with Eq. (jA.37p and us- 
ing in the latter the sum rule (jA.36p we obtain the 
condition 



IddiP,d) 



--Id^{a,l3,d) 



a=2f3-d+2 



(A.46) 



In d = 3 spatial dimensions this equation has the 
two solutions [13 



f3 



1, 



[3 ~ 0.795. 



(A.47) 



Only the first solution is found in angular approx- 
imation, see Eq. (jA.3ip . Both solutions are found 
numerically in the confined phase, see Refs. [l^ and 
(T9l |. respectively, 
ii) If x(^) is IR finite or vanishing the finite- 
temperature part of the tadpole [n] in the gap 
equation (|106p guarantees that uj{k) is IR finite but 
not vanishing, i.e. a = 0. In this case we find from 
the sum rule (|A.36p 



/3 



d-2 



(A.48) 



For d = 3 this yields /3 — ^, which is the solution 
realized in the deconfined phase, as our numerical 
solutions show. 
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